Sebastiaan Verbesselt
Instituut voor Natuur- en
Bosonderzoek
https://orcid.org/0000-0003-0173-1123
Tytti Jussila
Finnish Environment Institute
https://orcid.org/0000-0003-4646-0152
A datacube for the study area was already download from OpenEO based (see Retrieve Datacube OpenEO - Demervalley.Rmd script) for the period 01/01/2020 - 13/10/2025.
# Check automatically if you still need to install packages
list.of.packages <- c("tidyverse", "sf", "stars", "mapview", "lubridate", "dplyr", "rpart", "rpart.plot", "leaflet", "mapedit", "scales", "ggplot2", "rstudioapi","tidyr","zoo","np","kernlab","leafem","viridis","cubelyr","terra","signal","abind","INBOmd","INBOtheme")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
# Load the packages
library(tidyverse)
library(sf)
library(stars)
library(mapview)
library(lubridate)
library(dplyr)
library(rpart)
library(rpart.plot)
library(leaflet) # for interactive maps
library(leafem)
library(mapedit) # for drawing polygons interactively
library(scales)
library(ggplot2)
library(rstudioapi)
library(tidyr)
library(zoo)
library(np) # kernel regression
library(kernlab) # Gaussian processes
library(viridis)
library(terra)
library(signal)
library(abind)
library(INBOmd)
library(INBOtheme)
Refer to the location of your datafolder (interactively):
if (requireNamespace("rstudioapi", quietly = TRUE)) {
folder <- rstudioapi::selectDirectory(caption = "Select a folder")
print(folder)
} else {
message("rstudioapi not available; please install it with install.packages('rstudioapi').")
}
## NULL
if (is.null(folder)){
folder <- "./source/hydrology/data"
}
Load the datasets:
wetlands <- st_read(paste0(folder,"/raw/wetlands_kloosterbeemden.gpkg")) |> st_transform(32631) # Transform to local crs system (WGS 84 / UTM zone 31N - EPSG:32631)
## Reading layer `wetlands_kloosterbeemden' from data source
## `G:\Gedeelde drives\Team_BioDiv\5_Projecten\2024_Biodiversa_habitatpilot\WP2_3\INBO_Github_Biodiversa_Habitat_Pilot\Habitat_condition_indicators\source\hydrology\data\raw\wetlands_kloosterbeemden.gpkg'
## using driver `GPKG'
## Simple feature collection with 22 features and 43 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 191918.1 ymin: 188374.1 xmax: 193613.6 ymax: 189972.7
## Projected CRS: BD72 / Belgian Lambert 72
grasslands <- st_read(paste0(folder,"/raw/grasslands_kloosterbeemden.gpkg")) |> st_transform(32631) # Transform to local crs system (WGS 84 / UTM zone 31N - EPSG:32631)
## Reading layer `grasslands_kloosterbeemden' from data source
## `G:\Gedeelde drives\Team_BioDiv\5_Projecten\2024_Biodiversa_habitatpilot\WP2_3\INBO_Github_Biodiversa_Habitat_Pilot\Habitat_condition_indicators\source\hydrology\data\raw\grasslands_kloosterbeemden.gpkg'
## using driver `GPKG'
## Simple feature collection with 11 features and 43 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 192145.9 ymin: 188403 xmax: 193157.9 ymax: 189852.5
## Projected CRS: BD72 / Belgian Lambert 72
polygons <- rbind(wetlands,grasslands)
plot(polygons["HAB1"])
Select your area:
AOI <- polygons %>% dplyr::filter(HAB1 == "7140,rbbms")
plot(AOI["HAB1"])
# First: proxy loading (fast)
obj <- read_stars(paste0(folder,"/intermediate/polygon_Demervallei_1.nc"), proxy = TRUE) # Region code: 1: Kloosterbeemden, 2: Schulensmeer, 3: Webbekomsbroek
## B02, B03, B04, B05, B08, B8A, B11, B12, SCL,
obj2 <- read_stars(paste0(folder,"/intermediate/polygon_Demervallei_new_1.nc"), proxy = TRUE)
## B02, B03, B04, B05, B08, B8A, B11, B12, SCL,
obj3 <- read_stars(paste0(folder,"/intermediate/polygon_Demervallei_last_1.nc"), proxy = TRUE)
## B02, B03, B04, B05, B08, B8A, B11, B12, SCL,
# Ensure both have the same dimension names
# names(st_dimensions(obj))
# names(st_dimensions(obj2))
# names(st_dimensions(obj3))
# Drop any extra dimension (e.g., a band dimension named "X" or similar)
obj <- adrop(obj)
obj2 <- adrop(obj2)
obj3 <- adrop(obj3)
# Now combine along time
combined <- c(obj, obj2, along = "t")
combined <- c(combined, obj3, along = "t")
# st_dimensions(obj)
# st_dimensions(obj2)
# st_dimensions(obj3)
st_dimensions(combined)
## from to offset delta refsys values x/y
## x 1 383 636540 10 WGS 84 / UTM zone 31N NULL [x]
## y 1 282 5654200 -10 WGS 84 / UTM zone 31N NULL [y]
## t 1 858 NA NA POSIXct 2020-01-01,...,2025-10-13
Remove the data layers that we wont use to make sure we have enough free memory.
obj <- combined # Assign the combined datacube to the obj datacube (make a copy)
rm(obj2, obj3, combined) # Remove obj2, combined.
# Load now the data in memory:
obj <- st_as_stars(obj, along = "t")
wms_ortho_most_recent <- "https://geo.api.vlaanderen.be/OMWRGBMRVL/wms" # Most recent ortho image, winter, medium scale resolution.
wms_ortho <- "https://geo.api.vlaanderen.be/OMW/wms" # historical ortho images, winter, medium scale resolution.
wms_DEM <- "https://geo.api.vlaanderen.be/DHMV/wms" # Digital elevation model
wms_ANB <- "https://geo.api.vlaanderen.be/ANB/wms" # WMS layer ANB --> groenkaart (map of high vegetation)
Load the Jussila model
Jussila, Tytti, Risto K. Heikkinen, Saku Anttila, e.a. ‘Quantifying Wetness Variability in Aapa Mires with Sentinel-2: Towards Improved Monitoring of an EU Priority Habitat’. Remote Sensing in Ecology and Conservation 10, nr. 2 (2024): 172-87. https://doi.org/10.1002/rse2.363.
# Load the decision tree model
load(paste0(folder,"/raw/jussila_decisiontree.RData"))
# Visualize the decision tree structure
rpart.plot(tree_jussila, tweak = 1, extra = 0)
We will exclude “trees” from the image based on the most recent ortho image of Flanders. This because the inundation models are developed for “open” habitats. Normally, we could use regional datasets or European datasets for forests / tree cover / high vegetation cover to do the masking. E.g.: https://land.copernicus.eu/en/products/high-resolution-layer-forests-and-tree-cover. For Flanders, we masked trees interactively based on the High vegetation (> 3m) layer from the Flanders Groenkaart (2021) and the most recent available ortho image (winter, medium resolution).
outputfolder <- paste0(folder,"/processed")
This code (commented out) was used to create polygon shapes for the trees within the area. The final polygons are saved in an output directory.
# Ensure shapefile is in the right CRS (WGS84 lon/lat = EPSG:4326, since WCS expects that)
AOI_wgs84 <- st_transform(AOI, 4326)
# Extract bounding box
bb <- st_bbox(AOI_wgs84)
# Compute centroid of bbox
center_lng <- (bb["xmin"] + bb["xmax"]) / 2
center_lat <- (bb["ymin"] + bb["ymax"]) / 2
#
# map <- leaflet() %>%
# addWMSTiles(
# wms_ortho_most_recent,
# layers = "Ortho",
# options = WMSTileOptions(format = "image/png", transparent = FALSE), group = "Ortho") %>%
# addWMSTiles(
# wms_ANB,
# layers = "Grnkrt21",
# options = WMSTileOptions(format = "image/png", transparent = FALSE), group = "Grnkrt21") %>%
# fitBounds(
# lng1 = bb[["xmin"]],
# lat1 = bb[["ymin"]],
# lng2 = bb[["xmax"]],
# lat2 = bb[["ymax"]]
# ) %>%
# addPolygons(data = AOI_wgs84, color = "red", weight = 2) %>%
# addLayersControl(
# baseGroups = c("Ortho","Grnkrt21"),
# options = layersControlOptions(collapsed = FALSE)
# )
#
#
# # Allow interactive drawing of polygons and save as trees_object
# drawn <- mapedit::editMap(map)
#
# # This returns an sf object with drawn features
# trees_object <- drawn[["drawn"]]
# trees_object <- st_make_valid(trees_object)
# st_write(trees_object, paste0(outputfolder, "/", "trees_KB.shp"))
Now, we can mask out the trees form the study area.
# Inspect result
trees_object <- st_read(paste0(outputfolder, "/", "trees_KB.shp"))
## Reading layer `trees_KB' from data source
## `G:\Gedeelde drives\Team_BioDiv\5_Projecten\2024_Biodiversa_habitatpilot\WP2_3\INBO_Github_Biodiversa_Habitat_Pilot\Habitat_condition_indicators\source\hydrology\data\processed\trees_KB.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 7 features and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 4.986846 ymin: 51.01306 xmax: 4.990369 ymax: 51.01392
## Geodetic CRS: WGS 84
print(trees_object)
## Simple feature collection with 7 features and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 4.986846 ymin: 51.01306 xmax: 4.990369 ymax: 51.01392
## Geodetic CRS: WGS 84
## X_lflt_d ftr_typ geometry
## 1 318 polygon POLYGON ((4.987228 51.01332...
## 2 345 polygon POLYGON ((4.987841 51.01316...
## 3 379 polygon POLYGON ((4.988627 51.01315...
## 4 483 polygon POLYGON ((4.989897 51.01311...
## 5 549 polygon POLYGON ((4.989837 51.01349...
## 6 589 polygon POLYGON ((4.989322 51.01346...
## 7 613 polygon POLYGON ((4.98694 51.01349,...
par(mfrow=c(1,2))
plot(st_geometry(AOI_wgs84),border='grey',axes=T,main="Habitat boundary")
for (i in 1:nrow(trees_object)){
AOI_wgs84 <- st_difference(AOI_wgs84,trees_object[i,])
}
plot(st_geometry(AOI_wgs84,border='grey',axes=T,main="Habitat boundary without trees"))
par(mfrow=c(1,1))
See how inundation changes in in the past (based on winter ortho images) for your area. Select the year the want to visualize.
# Website with the public wms files for Flanders: https://wms.michelstuyts.be/?lang=nl
layer2024 <- "OMWRGB24VL"# select the winter image of 2024
layer2023 <- "OMWRGB23VL" # select the winter image of 2023
layer2022 <- "OMWRGB22VL"# select the winter image of 2022
layer2021 <- "OMWRGB21VL"# select the winter image of 2021
layer2020 <- "OMWRGB20VL"# select the winter image of 2020
map <- leaflet() %>%
addWMSTiles(
wms_ortho,
layers = layer2020,
options = WMSTileOptions(format = "image/png", transparent = FALSE),
group = "OMWRGB20VL"
) %>%
addWMSTiles(
wms_ortho,
layers = layer2021,
options = WMSTileOptions(format = "image/png", transparent = FALSE),
group = "OMWRGB21VL"
) %>%
addWMSTiles(
wms_ortho,
layers = layer2022,
options = WMSTileOptions(format = "image/png", transparent = FALSE),
group = "OMWRGB22VL"
) %>%
addWMSTiles(
wms_ortho,
layers = layer2023,
options = WMSTileOptions(format = "image/png", transparent = FALSE),
group = "OMWRGB23VL"
) %>%
addWMSTiles(
wms_ortho,
layers = layer2024,
options = WMSTileOptions(format = "image/png", transparent = FALSE),
group = "OMWRGB24VL"
) %>%
fitBounds(
lng1 = bb[["xmin"]],
lat1 = bb[["ymin"]],
lng2 = bb[["xmax"]],
lat2 = bb[["ymax"]]
) %>%
addPolygons(data = AOI_wgs84, color = "red", weight = 2) %>%
addLayersControl(
baseGroups = c("OMWRGB20VL","OMWRGB21VL","OMWRGB22VL","OMWRGB23VL","OMWRGB24V"),
options = layersControlOptions(collapsed = FALSE)
)
map
Most recent:
Date winter images study area:
most recent: 5/03/2025 (for Kloosterbeemden ), 28/03/2025 (for Schulensmeer), Webbekomsbroek (both 5/03/2025 (west side) and 28/03/2025 (east side).
2024: 06/04/2024 (for Kloosterbeemden & Webbekomsbroek), 30/07/2024 for Schulensmeer
2023: 01/03/2023 (for Kloosterbeemden & Webbekomsbroek), 31/05/2023 for Schulensmeer
2022: 04/03/2022 (for Kloosterbeemden ),18/03/2022 (for Schulensmeer), Webbekomsbroek (both 4/03/2025 (west side) and 18/03/2025 (east side).
2021: 01/03/2021 (for Kloosterbeemden & Webbekomsbroek), 21/02/2021 for Schulensmeer
2020: 24/03/2020 (for Kloosterbeemden & Webbekomsbroek), 17/03/2020 for Schulensmeer
Preprocess the datacube:
# Sentinel-2 SCL class codes and definitions
dplyr::tibble(
SCL = 0:11,
SCL_name = c(
"No data",
"Saturated or defective",
"Dark area pixels",
"Cloud shadow",
"Vegetation",
"Bare soils",
"Water",
"Cloud low probability / Unclassified",
"Cloud medium probability",
"Cloud high probability",
"Thin cirrus",
"Snow or ice"
)
) -> scl_code
print(scl_code)
## # A tibble: 12 × 2
## SCL SCL_name
## <int> <chr>
## 1 0 No data
## 2 1 Saturated or defective
## 3 2 Dark area pixels
## 4 3 Cloud shadow
## 5 4 Vegetation
## 6 5 Bare soils
## 7 6 Water
## 8 7 Cloud low probability / Unclassified
## 9 8 Cloud medium probability
## 10 9 Cloud high probability
## 11 10 Thin cirrus
## 12 11 Snow or ice
# Find all scenes with at least 95% of pixels classified as SCL 4, 5, or 6 (vegetation, bare soils, water)
# Mask out all remaining pixels with SCL values not equal to 4, 5, or 6
clear_dates <-
obj |>
as_tibble() |>
group_by(t) |>
summarize(prop_scl = sum(if_else(SCL %in% c(4,5,6), 1, 0)) / n()) |>
dplyr::filter(prop_scl > 0.80) |> # We can change this threshold if necessary.
pull(t)
obj_clear <-
obj |>
dplyr::filter(t %in% clear_dates) |>
mutate(across(everything(), ~ if_else(SCL %in% c(4, 5, 6), ., NA)))
rm(obj) # we won't use this object anymore.
# Re-project your area without trees back to the local Belgian crs system:
AOI <- st_transform(AOI_wgs84, 32631)
# Mask out the polygon
obj_poly <-
obj_clear |>
st_crop(AOI)
rm(obj_clear) # we won't use this object anymore.
# Convert stars object to a data frame for prediction
obj_df <- as.data.frame(obj_poly)
names(obj_df) <- c("x","y","t","b02","b03","b04","b05","b08","b8a","b11","b12","SCL")
# Add/calculate the necessary indices for the classification
obj_df$mndwi12 <- (obj_df$b03 - obj_df$b12) / (obj_df$b03 + obj_df$b12)
obj_df$mndwi11 <- (obj_df$b03 - obj_df$b11) / (obj_df$b03 + obj_df$b11)
obj_df$ndvi <- (obj_df$b8a - obj_df$b04) / (obj_df$b8a + obj_df$b04)
obj_df$ndwi_mf <- (obj_df$b03 - obj_df$b8a) / (obj_df$b03 + obj_df$b8a)
obj_df$ndmi_gao11 <- (obj_df$b8a - obj_df$b11) / (obj_df$b8a + obj_df$b11)
# additional indices. STR should be a good indication of moisture
swir_to_str <- function(swir) { # function to calculate moisture index STR (based on SWIR band 11 or 12)
swir <- swir/10000
STR <- ((1-swir)^2)/(2*swir) #5.29
return(STR)
}
obj_df$STR1 <- swir_to_str(obj_df$b11)
obj_df$STR2 <- swir_to_str(obj_df$b12)
summary(obj_df)
## x y t
## Min. :639375 Min. :5653165 Min. :2020-01-06 00:00:00
## 1st Qu.:639425 1st Qu.:5653185 1st Qu.:2021-02-28 06:00:00
## Median :639485 Median :5653210 Median :2022-07-26 12:00:00
## Mean :639485 Mean :5653210 Mean :2022-10-08 19:00:00
## 3rd Qu.:639545 3rd Qu.:5653235 3rd Qu.:2024-06-30 00:00:00
## Max. :639595 Max. :5653255 Max. :2025-10-10 00:00:00
##
## b02 b03 b04 b05
## Min. : -32.0 Min. : 52.0 Min. : 3.0 Min. : 176.0
## 1st Qu.: 256.0 1st Qu.: 417.0 1st Qu.: 311.0 1st Qu.: 742.0
## Median : 309.0 Median : 504.0 Median : 385.0 Median : 863.0
## Mean : 324.3 Mean : 502.2 Mean : 426.4 Mean : 865.8
## 3rd Qu.: 366.0 3rd Qu.: 583.0 3rd Qu.: 510.0 3rd Qu.: 991.0
## Max. :1846.0 Max. :1784.0 Max. :1576.0 Max. :1771.0
## NA's :21579 NA's :21579 NA's :21579 NA's :21579
## b08 b8a b11 b12
## Min. : 186 Min. : 152 Min. : 119 Min. : 56.0
## 1st Qu.:1667 1st Qu.:1804 1st Qu.:1440 1st Qu.: 727.0
## Median :2447 Median :2583 Median :1811 Median : 937.0
## Mean :2488 Mean :2593 Mean :1792 Mean : 993.3
## 3rd Qu.:3264 3rd Qu.:3396 3rd Qu.:2177 3rd Qu.:1201.0
## Max. :5669 Max. :5478 Max. :4105 Max. :2622.0
## NA's :21579 NA's :21579 NA's :21579 NA's :21579
## SCL mndwi12 mndwi11 ndvi
## Min. :4.000 Min. :-0.8365 Min. :-0.9100 Min. :-0.3214
## 1st Qu.:4.000 1st Qu.:-0.4187 1st Qu.:-0.6230 1st Qu.: 0.5726
## Median :4.000 Median :-0.3049 Median :-0.5581 Median : 0.7196
## Mean :4.211 Mean :-0.2996 Mean :-0.5389 Mean : 0.6783
## 3rd Qu.:4.000 3rd Qu.:-0.1914 3rd Qu.:-0.4882 3rd Qu.: 0.8150
## Max. :6.000 Max. : 0.6543 Max. : 0.4293 Max. : 0.9949
## NA's :21579 NA's :21579 NA's :21579 NA's :21579
## ndwi_mf ndmi_gao11 STR1 STR2
## Min. :-0.9374 Min. :-0.3420 Min. : 0.4233 Min. : 1.038
## 1st Qu.:-0.7280 1st Qu.: 0.0505 1st Qu.: 1.4056 1st Qu.: 3.223
## Median :-0.6756 Median : 0.1899 Median : 1.8515 Median : 4.383
## Mean :-0.6489 Mean : 0.1688 Mean : 2.5330 Mean : 5.258
## 3rd Qu.:-0.6029 3rd Qu.: 0.3127 3rd Qu.: 2.5442 3rd Qu.: 5.914
## Max. : 0.2475 Max. : 0.7513 Max. :41.0228 Max. :88.288
## NA's :21579 NA's :21579 NA's :21579 NA's :21579
Since the radiometric values of Sentinel-2 L2 products can have a different offset values (due to the change in processing baseline (4.00), introduces in January 2022): 0 or 1000. the DN (Digital Number) values have a radiometric offset of -1000 applied to ensure negative values can be represented. We check the min. and max. values within our datacube. https://sentiwiki.copernicus.eu/web/s2-processing#S2Processing-RadiometricOffset
# Check if there are no issues with the DN Values of the datacube (0-10000, so no offset with 1000).
# Check max values per column
(min_vals <- sapply(obj_df[,4:11], function(x) if(is.numeric(x)) min(x, na.rm = TRUE) else NA))
## b02 b03 b04 b05 b08 b8a b11 b12
## -32 52 3 176 186 152 119 56
(max_vals <- sapply(obj_df[,4:11], function(x) if(is.numeric(x)) max(x, na.rm = TRUE) else NA))
## b02 b03 b04 b05 b08 b8a b11 b12
## 1846 1784 1576 1771 5669 5478 4105 2622
# Warn if any numeric column has a minimum value that exceeds 1000
if (any(min_vals > 1000, na.rm = TRUE)) {
warning("Some columns in obj_df have min values greater than 1000")
}
# Warn if any numeric column exceeds 10000
if (any(max_vals > 10000, na.rm = TRUE)) {
warning("Some columns in obj_df have max values greater than 10000")
}
# Keep the original x, y, t dimension form your stars object.
# If we convert our stars object to a dataframe for classification, we lose the x, y, t information. When we want to convert it back to an array/datacube, we need to know the original dimensions.
(dim <- st_dimensions(obj_poly)) # See the dimensions of x y and t
## from to offset delta refsys values x/y
## x 284 306 636540 10 WGS 84 / UTM zone 31N NULL [x]
## y 95 104 5654200 -10 WGS 84 / UTM zone 31N NULL [y]
## t 1 168 NA NA POSIXct 2020-01-06,...,2025-10-10
nx <- dim$x$to - dim$x$from + 1 # size in the x dimension
ny <- dim$y$to - dim$y$from + 1 # size in the y dimension
nt <- dim$t$to - dim$t$from + 1 # size in the t dimension
predictions <- predict(tree_jussila, obj_df, type = "class")
# Convert factor to numeric codes
pred_numeric <- as.numeric(predictions)
# Set predictions to NA if the corresponding row in obj_df contains any NA (masked out values are set to NA instead of 1 or 2)
pred_numeric[!complete.cases(obj_df)] <- NA
# Store levels for later labeling
(pred_levels <- levels(predictions))
## [1] "dry" "water"
# Reshape to array
pred_array <- array(pred_numeric, dim = c(nx, ny, nt))
# Create stars object
obj_classified <- st_as_stars(pred_array, dimensions = st_dimensions(obj_poly))
st_crs(obj_classified) <- st_crs(obj_poly)
# Visualize
# Convert the stars object to factor using st_set_dimensions
obj_classified_factor <- obj_classified[AOI]
obj_classified_factor[[1]] <- factor(
obj_classified_factor[[1]],
levels = c(1, 2),
labels = pred_levels
)
# Plot with discrete scale
ggplot() +
ggtitle("Class. (Jussila)") +
geom_stars(data = obj_classified_factor) +
geom_sf(data = AOI, fill = NA, color = "red", linewidth = 1) +
facet_wrap(~t, ncol=25) +
theme_minimal() +
theme(
strip.text = element_text(size = 6), # <-- make facet titles smaller,
axis.text = element_blank(), # remove axis tick labels
axis.ticks = element_blank(), # remove tick marks
) +
scale_fill_manual(
name = "Class",
values = c("dry" = "tan", "water" = "blue")
)
Lefebvre, Gaëtan, Aurélie Davranche, Loïc Willm, e.a. ‘Introducing WIW for Detecting the Presence of Water in Wetlands with Landsat and Sentinel Satellites’. Remote Sensing 11, nr. 19 (2019): 2210. https://doi.org/10.3390/rs11192210.
watercl_wiw <- obj_df$b8a <= 1804 & obj_df$b12 <= 1131
watercl_wiw <- watercl_wiw*1
pred_levels <- c("dry" , "water")
# Reshape to array
wiw_array <- array(watercl_wiw, dim = c(nx, ny, nt))
# Create stars object
obj_classified_wiw <- st_as_stars(wiw_array, dimensions = st_dimensions(obj_poly))
st_crs(obj_classified_wiw) <- st_crs(obj_poly)
# Visualize
# Convert the stars object to factor using st_set_dimensions
obj_classified_wiw_factor <- obj_classified_wiw[AOI]
obj_classified_wiw_factor[[1]] <- factor(
obj_classified_wiw_factor[[1]],
levels = c(0, 1),
labels = pred_levels
)
# Plot with discrete scale
ggplot() +
ggtitle("Class. (WiW)") +
geom_stars(data = obj_classified_wiw_factor) +
geom_sf(data = AOI, fill = NA, color = "red", linewidth = 1) +
facet_wrap(~t, ncol=25) +
theme_minimal()+
theme(
strip.text = element_text(size = 6), # <-- make facet titles smaller,
axis.text = element_blank(), # remove axis tick labels
axis.ticks = element_blank(), # remove tick marks
) +
scale_fill_manual(
name = "Class",
values = c("dry" = "tan", "water" = "blue")
)
We can check the output of both classifications and compare it with an ortho image with the same date. E.g. 01/03/2023.
# Extract the time values
times <- st_get_dimension_values(obj_classified_wiw, "t")
# Define your desired time
target_time <- as.POSIXct("2023-03-01", tz = "UTC")
# Find index of the nearest timestamp
i <- which.min(abs(times - target_time))
# Extract that layer
layer_2023_03_01_wiw <- obj_classified_wiw[,,,i]
layer_2023_03_01_Tytti <- obj_classified[,,,i] - 1 # convert values [1,2] --> [0,1]
layer_2023_03_01_wiw <- layer_2023_03_01_wiw %>% st_warp(crs=4326)
layer_2023_03_01_Tytti <- layer_2023_03_01_Tytti %>% st_warp(crs=4326)
# Define your color palette
palette_function <- function(stars_data){
if (min(stars_data, na.rm=T) == max(stars_data, na.rm=T)){
if (max(stars_data, na.rm=T)==1){
palette <- colorFactor(palette = "blue",domain = 1,na.color = "transparent")
} else {
palette <- colorFactor(palette = "tan",domain = 0,na.color = "transparent")
}
}else {
palette <- colorFactor(palette = c("tan","blue"),domain= c(0,1),na.color = "transparent")
}
}
# Base map (your existing code)
layer <- "OMWRGB23VL"
map <- leaflet() %>%
addTiles(group = "OpenStreetMap") %>% # <-- Adds OSM as base map
addWMSTiles(
wms_ortho,
layers = layer,
options = WMSTileOptions(format = "image/png", transparent = FALSE),
group = "OMWRGB20VL"
) %>%
addWMSTiles(
wms_DEM,
layers = "DHMV_II_HILL_25cm",
options = WMSTileOptions(format = "image/png", transparent = TRUE),
group = "DHMV_II_HILL_25cm"
) %>%
fitBounds(
lng1 = bb[["xmin"]],
lat1 = bb[["ymin"]],
lng2 = bb[["xmax"]],
lat2 = bb[["ymax"]]
) %>%
addPolygons(data = AOI_wgs84, color = "red", weight = 2) %>%
addStarsImage(
layer_2023_03_01_wiw,
colors = palette_function(layer_2023_03_01_wiw[["X"]]),
opacity = 0.8,
group = "2023-03-01 Wiw"
) %>%
addStarsImage(
layer_2023_03_01_Tytti,
colors = palette_function(layer_2023_03_01_Tytti[["X"]]),
opacity = 0.8,
group = "2023-03-01 Tytti"
) %>%
addLayersControl(
baseGroups = c("OpenStreetMap","OMWRGB23VL", "DHMV_II_HILL_25cm"),
overlayGroups = c("2023-03-01 Wiw", "2023-03-01 Tytti"),
options = layersControlOptions(collapsed = FALSE)
)
map
Add soil moisture/wetness indices to your stars object.
obj_wetness <-
obj_poly |>
mutate(NDMI = (B8A - B11) / (B8A + B11), # Gao's Moisture index
NDWI = (B03 - B8A) / (B03 + B8A), # Mcfeeters Water index
NDPI = (B03-B11)/(B03+B11), # Pond index. from slovakia Clizsky potok example
STR = ((1-B11/10000)^2)/(2*B11/10000)) # Transformed SWIR. Should be linearly correlated with soil moisture (Sadeghi et al.,2017, https://doi.org/10.1016/j.rse.2017.05.041)
Visualize the different soil moisture indices for the target date.
# Extract that layer
obj_wetness_target <- obj_wetness[,,,i]
obj_wetness_target <- obj_wetness_target %>% st_warp(crs=4326)
# Extract the time value from the layer
date_value <- st_get_dimension_values(obj_wetness_target, "t")
date_value
## [1] "2023-03-01 UTC"
# Define a viridis palette (yellow -> blue)
# direction = -1 flips the default viridis (blue -> yellow) to yellow -> blue
viridis_pal <- viridis::viridis(256, option = "D", direction = -1)
# You can define a function to generate color scales for each layer dynamically
palette_function <- function(layer) {
# Extract the numeric values safely
layer_values <- as.vector(st_as_stars(layer)[[1]])
leaflet::colorNumeric(
palette = viridis_pal,
domain = layer_values,
na.color = "transparent"
)
}
# Create map with OSM background
map <- leaflet() %>%
addTiles(group = "OpenStreetMap") %>%
addWMSTiles(
wms_ortho,
layers = layer,
options = WMSTileOptions(format = "image/png", transparent = FALSE),
group = "OMWRGB20VL"
) %>%
addWMSTiles(
wms_DEM,
layers = "DHMV_II_HILL_25cm",
options = WMSTileOptions(format = "image/png", transparent = TRUE),
group = "DHMV_II_HILL_25cm"
) %>%
fitBounds(
lng1 = bb[["xmin"]],
lat1 = bb[["ymin"]],
lng2 = bb[["xmax"]],
lat2 = bb[["ymax"]]
) %>%
addPolygons(data = AOI_wgs84, color = "red", weight = 2) %>%
addStarsImage(
-obj_wetness_target["B11"],
colors = palette_function(-obj_wetness_target["B11"]), # multiply with -1 (because high values normally inidicate dry spots, not wet spots like in the other indices).
opacity = 0.8,
group = paste0(date_value, " B11")
) %>%
addStarsImage(
obj_wetness_target["NDWI"],
colors = palette_function(obj_wetness_target["NDWI"]),
opacity = 0.8,
group = paste0(date_value, " NDWI")
) %>%
addStarsImage(
obj_wetness_target["NDMI"],
colors = palette_function(obj_wetness_target["NDMI"]),
opacity = 0.8,
group = paste0(date_value, " NDMI")
) %>%
addStarsImage(
obj_wetness_target["NDPI"],
colors = palette_function(obj_wetness_target["NDPI"]),
opacity = 0.8,
group = paste0(date_value, " NDPI")
) %>%
addStarsImage(
obj_wetness_target["STR"],
colors = palette_function(obj_wetness_target["STR"]),
opacity = 0.8,
group = paste0(date_value, " STR")
) %>%
addLayersControl(
baseGroups = c("OpenStreetMap","OMWRGB20VL", "DHMV_II_HILL_25cm"),
overlayGroups = c(
paste0(date_value, " B11"),
paste0(date_value, " NDWI"),
paste0(date_value, " NDMI"),
paste0(date_value, " NDPI"),
paste0(date_value, " STR")
),
options = layersControlOptions(collapsed = FALSE)
)
map
Juss_cl <- obj_classified-1 # convert values form 1-2 to 0-1
obj_wetness <- c(obj_wetness,Juss_cl)
obj_wetness <- c(obj_wetness,obj_classified_factor)
obj_wetness <- c(obj_wetness,obj_classified_wiw)
obj_wetness <- c(obj_wetness,obj_classified_wiw_factor)
names(obj_wetness)[14] <- "Juss_cl"
names(obj_wetness)[15] <- "Juss_cl_fact"
names(obj_wetness)[16] <- "Wiw_cl"
names(obj_wetness)[17] <- "Wiw_cl_fact"
obj_wetness
## stars object with 3 dimensions and 17 attributes
## attribute(s):
## B02 B03 B04 B05
## Min. : -32.0 Min. : 52.0 Min. : 3.0 Min. : 176.0
## 1st Qu.: 256.0 1st Qu.: 417.0 1st Qu.: 311.0 1st Qu.: 742.0
## Median : 309.0 Median : 504.0 Median : 385.0 Median : 863.0
## Mean : 324.3 Mean : 502.2 Mean : 426.4 Mean : 865.8
## 3rd Qu.: 366.0 3rd Qu.: 583.0 3rd Qu.: 510.0 3rd Qu.: 991.0
## Max. :1846.0 Max. :1784.0 Max. :1576.0 Max. :1771.0
## NA's :21579 NA's :21579 NA's :21579 NA's :21579
## B08 B8A B11 B12
## Min. : 186 Min. : 152 Min. : 119 Min. : 56.0
## 1st Qu.:1667 1st Qu.:1804 1st Qu.:1440 1st Qu.: 727.0
## Median :2447 Median :2583 Median :1811 Median : 937.0
## Mean :2488 Mean :2593 Mean :1792 Mean : 993.3
## 3rd Qu.:3264 3rd Qu.:3396 3rd Qu.:2177 3rd Qu.:1201.0
## Max. :5669 Max. :5478 Max. :4105 Max. :2622.0
## NA's :21579 NA's :21579 NA's :21579 NA's :21579
## SCL NDMI NDWI NDPI
## Min. :4.000 Min. :-0.3420 Min. :-0.9374 Min. :-0.9100
## 1st Qu.:4.000 1st Qu.: 0.0505 1st Qu.:-0.7280 1st Qu.:-0.6230
## Median :4.000 Median : 0.1899 Median :-0.6756 Median :-0.5581
## Mean :4.211 Mean : 0.1688 Mean :-0.6489 Mean :-0.5389
## 3rd Qu.:4.000 3rd Qu.: 0.3127 3rd Qu.:-0.6029 3rd Qu.:-0.4882
## Max. :6.000 Max. : 0.7513 Max. : 0.2475 Max. : 0.4293
## NA's :21579 NA's :21579 NA's :21579 NA's :21579
## STR Juss_cl Juss_cl_fact Wiw_cl Wiw_cl_fact
## Min. : 0.4233 Min. :0.000 dry :13154 Min. :0.0000 dry :13950
## 1st Qu.: 1.4056 1st Qu.:0.000 water: 3907 1st Qu.:0.0000 water: 3111
## Median : 1.8515 Median :0.000 NA's :21579 Median :0.0000 NA's :21579
## Mean : 2.5330 Mean :0.229 Mean :0.1823
## 3rd Qu.: 2.5442 3rd Qu.:0.000 3rd Qu.:0.0000
## Max. :41.0228 Max. :1.000 Max. :1.0000
## NA's :21579 NA's :21579 NA's :21579
## dimension(s):
## from to offset delta refsys values x/y
## x 284 306 636540 10 WGS 84 / UTM zone 31N NULL [x]
## y 95 104 5654200 -10 WGS 84 / UTM zone 31N NULL [y]
## t 1 168 NA NA POSIXct 2020-01-06,...,2025-10-10
## Temporal aggregation: 2weekly median value with customized approach
# generate 2-week breaks for the observed years
time_vals <- st_get_dimension_values(obj_wetness, "t")
start <- floor_date(min(time_vals), "month")
end <- ceiling_date(max(time_vals), "month")
# All 1st and 15th of each month between start and end
breaks <- sort(c(seq(start, end, by = "1 month"),
seq(start + days(14), end, by = "1 month")))
# Aggregate using 2-week intervals
obj_2week <- aggregate(obj_wetness, by = breaks, FUN = median, na.rm = TRUE)
# Summarise 2-weekly
table_2w <-
obj_2week[,,,] |>
as_tibble() |>
group_by(time) |>
summarize(inundation = mean(Juss_cl, na.rm = TRUE)*100,
inundation_Wiw = mean(Wiw_cl, na.rm = TRUE)*100,
NDMI = mean(NDMI, na.rm=T), #Ranges from -1 to +1, where the lowest values indicate low vegetation water content, and the highest ones correspond to high water content
NDWI = mean(NDWI, na.rm=T), #Over 0.2: water. Below 0: non-water
NDPI = mean(NDPI, na.rm=T), # scores above a certain high value, e.g., larger than 0.75, indicates a pond area with good water quality
STR = mean(STR, na.rm=T), # Higher values, higher soil moisture content
B11 = mean(B11, na.rm=T)) # lower values, higher moisture.
table_2w_plotting <- table_2w[!is.na(table_2w$inundation),] # drop NAs for plotting
ggplot(table_2w_plotting) + ## check the dates in table to plot a specific single year
aes(x = time, y = inundation) +
ylim(0,100) +
geom_line() + geom_point() +
labs(title = "Inundated area %",
x = "Date",
y = "Inundated area (%)") +
theme_minimal()
ggplot(table_2w_plotting) +
aes(x = time, y = STR) +
#ylim(1,5) +
geom_line() + geom_point() +
labs(title = "STR",
x = "Date",
y = "STR") +
theme_minimal()
The time series is still irregular. We will use linear interpolation for gap filling.
# extract 2-weekly time points
time_vals2w <- st_get_dimension_values(obj_2week, "time")
# interpolate NA values along time dimension
obj_filled_2w <- st_apply(
obj_2week,
MARGIN = c("x", "y"),
FUN = function(ts) { # repeat interpolation function over each pixel time series
if (all(is.na(ts))) { # keep NA time series as NA outside site polygon
return(ts)
}
approx( # interpolate missing time points
x = as.numeric(time_vals2w),
y = ts,
xout = as.numeric(time_vals2w),
method = "linear",
rule = 2
)$y
},.fname = "time"
)
# fix the broken time dimension in output
obj_filled_2w <- st_set_dimensions(obj_filled_2w, "time", values = time_vals2w)
ggplot() +
ggtitle("2-weekly mean class (Jussila)") +
geom_stars(data = obj_filled_2w["Juss_cl"]) +
geom_sf(data = AOI, fill = NA, color = "red", linewidth = 1) +
facet_wrap(~time,ncol=20) +
theme_minimal() +
theme(axis.text = element_blank(), # remove axis tick labels
axis.ticks = element_blank(), # remove tick marks
) +
scale_fill_gradientn(
name = "Wetness",
colors = c("tan", "cyan", "blue"),
na.value = "gray",
limits = c(0, 1)
)
# Plot 2-week time series (gapfilled)
# first summarise for site area:
table_gf2w <- # Gapfilled 2-weekly table
obj_filled_2w[,,,] |>
as_tibble() |>
group_by(time) |>
summarize(inundation = mean(Juss_cl, na.rm = TRUE)*100,
inundation_Wiw = mean(Wiw_cl, na.rm = TRUE)*100,
NDMI = mean(NDMI, na.rm=T), #Ranges from -1 to +1, where the lowest values indicate low vegetation water content, and the highest ones correspond to high water content
NDWI = mean(NDWI, na.rm=T), #Over 0.2: water. Below 0: non-water
NDPI = mean(NDPI, na.rm=T), # scores above a certain high value, e.g., larger than 0.75, indicates a pond area with good water quality
STR = mean(STR, na.rm=T), # Higher values, higher soil moisture content
B11 = mean(B11, na.rm=T)) # lower values, higher moisture.
# New column: Flag gapfilled values
table_gf2w$gapfilled <- "2w median"
table_gf2w$gapfilled[is.na(table_2w$B11)] <- "gapfilled"
table_gf2w$gapfilled <- as.factor(table_gf2w$gapfilled)
# time format to Date for plotting
table_gf2w$date <- as.Date(table_gf2w$time)
ggplot(table_gf2w) + #check the dates in table to plot a single year ([65:74,] -> 2018 in this case)
aes(x = date, y = inundation) + scale_x_date(date_labels = "%b %y") +
geom_line() + geom_point(aes(colour = gapfilled), size=2) +
scale_color_manual(values = c("gapfilled" = "red", "2w median" = "black"), name=NULL)+
labs(title = "Inundation %",x = "Date",
y = "Inundated area (%)") +
ylim(0,100)+
theme_minimal()
ggplot(table_gf2w) +
aes(x = date, y = STR) + scale_x_date(date_labels = "%b %y") +
geom_line() + geom_point(aes(colour = gapfilled), size=2) +
scale_color_manual(values = c("gapfilled" = "red", "2w median" = "black"), name=NULL)+
labs(title = "STR", x = "Date",
y = "STR") +
#ylim(1,5)+ # adjust if needed
theme_minimal()
ggplot(table_gf2w) +
aes(x = date, y = NDMI) + scale_x_date(date_labels = "%b %y") +
geom_line() + geom_point(aes(colour = gapfilled), size=2) +
scale_color_manual(values = c("gapfilled" = "red", "2w median" = "black"), name=NULL)+
labs(title = "NDMI", x = "Date",
y = "NDMI") +
ylim(-0.1, 0.55)+ # adjust if needed
theme_minimal()
ggplot(table_gf2w) +
aes(x = date, y = -B11) + scale_x_date(date_labels = "%b %y") +
geom_line() + geom_point(aes(colour = gapfilled), size=2) +
scale_color_manual(values = c("gapfilled" = "red", "2w median" = "black"), name=NULL)+
labs(title = "SWIR",x = "Date",
y = "B11 (neg)") + # plot SWIR as negation for comparability (low SWIR, high moisture. Opposite to other indicators)
#ylim(-2500, -800)+ # adjust if needed
theme_minimal()
Plotting yearly indicators.
# Using 2-weekly gapfilled stars object to calculate statistics
# aggregate stars into yearly layers
## Mean (for each pixel?)
obj_y_mean <- aggregate(obj_filled_2w, by = "year",
FUN=mean, na.rm=T)
## Min (for each pixel?)
obj_y_min <- aggregate(obj_filled_2w, by = "year",
FUN=min, na.rm=T) # returns Inf for NA areas
obj_y_min[obj_y_min == Inf] <- NA # convert Inf back to NA
## Max (for each pixel?)
obj_y_max <- aggregate(obj_filled_2w, by = "year",
FUN=max, na.rm=T) # returns -Inf for NA areas
obj_y_max[obj_y_max == -Inf] <- NA # convert -Inf back to NA
## Plot yearly raster maps
# SWIR [7], inundation Jussila [14], inundation Wiw [16], NDMI = [10], NDWI = [11], NDPI = [12], STR = [13]
# Inundation (Jussila). Inundated % of time
plot(obj_y_mean[14], col = viridis(100, option = "D", direction=-1),
breaks = seq(0, 1, length.out = 101), main = "Relative frequency of inundated pixels by Jussila model")
plot(obj_y_mean[14]*26, col = viridis(100, option = "D", direction=-1),
breaks = seq(0, 26, length.out = 101), main = "Number of 2-weekly inundation per year by Jussila model")
# Inundation (Wiw). Inundated % of time
plot(obj_y_mean[16], col = viridis(100, option = "D", direction=-1),
breaks = seq(0, 1, length.out = 101), main = "Relative frequency of inundated pixels by Wiw model")
plot(obj_y_mean[16]*26, col = viridis(100, option = "D", direction=-1),
breaks = seq(0, 26, length.out = 101), main = "Number of 2-weekly inundation per year by Wiw model")
# STR
plot(obj_y_mean[13], col = viridis(100, option = "D", direction=-1),
breaks = seq(0, 40, length.out = 101))
# NDMI
plot(obj_y_mean[10], col = viridis(100, option = "D", direction= -1),
breaks = seq(min(obj_y_mean[[10]], na.rm=T),
max(obj_y_mean[[10]], na.rm=T), length.out = 101))
# Extract the layer
x <- obj_y_mean[14]
# Classify pixel values
vals <- x[[1]]
classified_vals <- ifelse(
is.na(vals), NA,
ifelse(vals == 0, "permanently dry",
ifelse(vals == 1, "permanently inundated",
"seasonally inundated"))
)
# Rebuild a stars object, preserving geometry
classified <- st_as_stars(classified_vals)
st_dimensions(classified) <- st_dimensions(x)
st_crs(classified) <- st_crs(x)
# Convert to factor (only with existing levels)
present_classes <- sort(unique(na.omit(as.vector(classified_vals))))
classified[[1]] <- factor(classified[[1]], levels = present_classes)
# Define full color scheme, but only use what’s needed
full_colors <- c(
"permanently dry" = "tan",
"seasonally inundated" = "lightblue",
"permanently inundated" = "blue4"
)
class_colors <- full_colors[present_classes]
ggplot() +
geom_stars(data = classified, aes(x = x, y = y)) +
scale_fill_manual(values = class_colors, name = "Inundation Class") +
geom_sf(data = AOI, fill = NA, color = "red", size = 0.7) +
coord_sf() +
facet_wrap(~ time, ncol = 3) +
labs(
title = "Inundation Frequency Classification",
x = "Longitude",
y = "Latitude"
) +
theme_minimal() +
theme(
legend.position = "right",
panel.border = element_rect(color = "grey40", fill = NA),
strip.text = element_text(face = "bold")
)
Plotting yearly indicators for site area:
# MEAN - summarise values for site area
table_y_mean <-
obj_y_mean[,,,] |>
as_tibble() |>
group_by(time) |>
summarize(inundation = mean(Juss_cl, na.rm = TRUE)*100,
inundation_Wiw = mean(Wiw_cl, na.rm = TRUE)*100,
NDMI = mean(NDMI, na.rm=T), #Ranges from -1 to +1, where the lowest values indicate low vegetation water content, and the highest ones correspond to high water content
NDWI = mean(NDWI, na.rm=T), #Over 0.2: water. Below 0: non-water
NDPI = mean(NDPI, na.rm=T), # scores above a certain high value, e.g., larger than 0.75, indicates a pond area with good water quality
STR = mean(STR, na.rm=T), # Higher values, higher soil moisture content
B11 = mean(B11, na.rm=T)) # lower values, higher moisture.
#summarise min for sites (site mean at minimum wet situation)
table_y_min <-
obj_y_min[,,,] |>
as_tibble() |>
group_by(time) |>
summarize(inundation = mean(Juss_cl, na.rm = TRUE)*100,
inundation_Wiw = mean(Wiw_cl, na.rm = TRUE)*100,
NDMI = mean(NDMI, na.rm=T), #Ranges from -1 to +1, where the lowest values indicate low vegetation water content, and the highest ones correspond to high water content
NDWI = mean(NDWI, na.rm=T), #Over 0.2: water. Below 0: non-water
NDPI = mean(NDPI, na.rm=T), # scores above a certain high value, e.g., larger than 0.75, indicates a pond area with good water quality
STR = mean(STR, na.rm=T), # Higher values, higher soil moisture content
B11 = mean(B11, na.rm=T)) # lower values, higher moisture.
# for SWIR, wetness minimum is SWIR max value.
#summarise max for sites (site mean at minimum wet situation)
table_y_max <-
obj_y_max[,,,] |>
as_tibble() |>
group_by(time) |>
summarize(inundation = mean(Juss_cl, na.rm = TRUE)*100,
inundation_Wiw = mean(Wiw_cl, na.rm = TRUE)*100,
NDMI = mean(NDMI, na.rm=T), #Ranges from -1 to +1, where the lowest values indicate low vegetation water content, and the highest ones correspond to high water content
NDWI = mean(NDWI, na.rm=T), #Over 0.2: water. Below 0: non-water
NDPI = mean(NDPI, na.rm=T), # scores above a certain high value, e.g., larger than 0.75, indicates a pond area with good water quality
STR = mean(STR, na.rm=T), # Higher values, higher soil moisture content
B11 = mean(B11, na.rm=T)) # lower values, higher moisture.
# for SWIR, wetness maximum is SWIR min value.
# plot all yearly stats in one graph: Mean and shaded range between min-max
# STR
ggplot() +
#ylim(1,5)+ # adjust if needed
#geom_line(data= table_gf2w, aes(x = time - months(6), y = STR), color="grey") + # full time series to background?
geom_line(data= table_y_mean, aes(x = time, y = STR)) +
geom_point(data= table_y_mean, aes(x = time, y = STR)) +
geom_ribbon(aes(x = table_y_mean$time,
ymin = table_y_min$STR,
ymax = table_y_max$STR),fill="#1f9e89", alpha= 0.2)+
labs(title = "STR moisture mean and range (min, max)",
x = "Year",
y = "STR") +
theme_minimal() + theme(legend.position="none")
# NDMI
ggplot() +
#ylim(-0.1,0.5)+ # adjust if needed
#geom_line(data= table_gf2w, aes(x = time - months(6), y = NDMI), color="grey") + # full time series to background?
geom_line(data= table_y_mean, aes(x = time, y = NDMI)) +
geom_point(data= table_y_mean, aes(x = time, y = NDMI)) +
geom_ribbon(aes(x = table_y_mean$time,
ymin = table_y_min$NDMI,
ymax = table_y_max$NDMI),fill="#1f9e89", alpha= 0.2)+
labs(title = "NDMI moisture mean and range (min, max)",
x = "Year",
y = "NDMI") +
theme_minimal() + theme(legend.position="none")
# plot percentage of permanently and seasonally wet area.
# Jussila model
ggplot() +
ylim(0,100)+
geom_ribbon(aes(x = table_y_mean$time, # permanent
ymin = 0,
ymax = table_y_min$inundation),fill="#1f9e89", alpha= 0.6)+
geom_ribbon(aes(x = table_y_mean$time, # seasonal
ymin = table_y_min$inundation,
ymax = table_y_max$inundation),fill="#6ece58", alpha= 0.3)+ #Temporarily inundated
geom_ribbon(aes(x = table_y_mean$time, # never
ymin = table_y_max$inundation,
ymax = 100),fill="#fde725", alpha= 0.2)+ # Never inundated
labs(title = "Permanently and seasonally inundated area (Jussila model)",
x = "Year",
y = "Inundated area %") +
theme_minimal() + theme(legend.position="none") # + scale_x_date(date_breaks = "1 year")
# Lefebvre Wiw model
ggplot() +
ylim(0,100)+
geom_ribbon(aes(x = table_y_mean$time, # permanent
ymin = 0,
ymax = table_y_min$inundation_Wiw),fill="#1f9e89", alpha= 0.6)+
geom_ribbon(aes(x = table_y_mean$time, # seasonal
ymin = table_y_min$inundation_Wiw,
ymax = table_y_max$inundation_Wiw),fill="#6ece58", alpha= 0.3)+
geom_ribbon(aes(x = table_y_mean$time, # never
ymin = table_y_max$inundation_Wiw,
ymax = 100),fill="#fde725", alpha= 0.2)+
labs(title = "Permanently and seasonally inundated area (Lefebvre Wiw)",
x = "Year",
y = "Inundated area %") +
theme_minimal() + theme(legend.position="none") # + scale_x_date(date_breaks = "1 year")
# Assuming your object is named obj_classified and has a time dimension "t"
st <- obj_classified[AOI]
# 1. Extract time dimension
t_existing <- st_get_dimension_values(st, "t")
# 2. Create complete 5-day sequence
t_full <- seq(min(t_existing), max(t_existing), by = "5 days")
# 3. Find missing dates
t_missing <- setdiff(t_full, t_existing)
# 4. Proceed only if there are missing timestamps
if (length(t_missing) > 0) {
# --- Template Creation: Clean Reset of Time Dimension ---
# Take the first time slice, retaining all 4 dimensions (x, y, band, t)
template <- st[,,,1, drop = FALSE]
# Replace data with NA
template[[1]][] <- NA
# Define a dummy date placeholder, ensuring it's a POSIXct object
dummy_date <- min(t_existing) - as.difftime(1, units = "days")
dummy_date <- as.POSIXct(dummy_date)
#Reset the time dimension of the template using st_set_dimensions.
# This overwrites the existing dimension value and ensures the internal structure
# of the dimension object remains valid (e.g., handles "intervals" properties).
template <- st_set_dimensions(
template,
"t",
values = dummy_date,
# Also reset the point flag to FALSE for consistency if using intervals
point = FALSE
)
# --- Creating and Assigning Missing Layers ---
# Create a list of new layers for missing timestamps
missing_layers <- lapply(t_missing, function(tt) {
new_layer <- template
# Set the correct missing date 'tt' (overwriting the dummy_date)
# The crucial fix from previous steps: assign the result back
new_layer <- st_set_dimensions(new_layer, "t", values = tt)
return(new_layer)
})
# --- Combining Layers ---
# Prepare the full list of objects to combine
all_layers_to_combine <- c(list(st), missing_layers)
# Combine all layers along the time dimension
# This uses the c.stars method correctly by unpacking the list
st_full <- do.call(c, c(all_layers_to_combine, along = "t"))
# --- Final Sort (To Interleave the NA layers and resolve stacking) ---
# Sort the time dimension chronologically.
t_values_full <- st_get_dimension_values(st_full, "t")
sorted_indices <- order(t_values_full)
# Reorder the object data by subsetting with the sorted time indices
st_full <- st_full[,,, sorted_indices, drop = FALSE]
# Ensure the dimension values are also updated to the sorted order
st_full <- st_set_dimensions(
st_full,
"t",
values = sort(t_values_full)
)
} else {
st_full <- st
}
# Result
st_full
## stars object with 3 dimensions and 1 attribute
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## X 1 1 1 1.229002 1 2 99319
## dimension(s):
## from to offset delta refsys values x/y
## x 284 306 636540 10 WGS 84 / UTM zone 31N NULL [x]
## y 95 104 5654200 -10 WGS 84 / UTM zone 31N NULL [y]
## t 1 506 NA NA POSIXct 2020-01-06,...,2025-10-10
st_ones_alt <- st_full
st_ones_alt[[1]][!is.na(st_ones_alt[[1]])] <- 1
st_ones_alt[[1]][is.na(st_ones_alt[[1]])] <- NA # Ensure NA remains NA (as the first method does)
# Aggregate st_full by month, calculating the mean over each month.
# FUN = mean calculates the sum and divides by the number of non-NA dates in the month.
st_monthly_sum <- aggregate(
st_ones_alt,
by = "months",
FUN = sum,
na.rm = TRUE
)
ggplot() +
ggtitle("sum of valid pixels per month") +
geom_stars(data = st_monthly_sum) +
geom_sf(data = AOI, fill = NA, color = "red", linewidth = 1) +
facet_wrap(~time,ncol=12) +
theme_minimal() +
theme(axis.text = element_blank(), # remove axis tick labels
axis.ticks = element_blank(), # remove tick marks
)+
scale_fill_gradientn(
name = "X",
colors = c("white", "purple","darkblue"),
na.value = "gray",
#limits = c(0, 1)
)
# --- STEP 1: Extract Spatial Template ---
# The function needs the extent and CRS to turn the plain array data back into a SpatRaster.
# Slice the first time step and convert it to a SpatRaster to get the template metadata.
# template raster for extent and CRS
r_template <- rast(slice(obj_classified, "t", 1))
template_crs <- crs(r_template, proj=TRUE)
template_extent <- ext(r_template)
# function to process a single slice
focal_na_fill_modal <- function(x_slice) {
r_slice <- rast(t(x_slice)) # transpose first
ext(r_slice) <- template_extent
crs(r_slice) <- template_crs
r_focal <- focal(
r_slice,
w = 3,
fun = modal,
na.policy = "only",
na.rm = TRUE
)
as.matrix(r_focal) |> t() # transpose back
}
# Apply across time
obj_focal_filled <- st_apply(obj_classified, MARGIN = "t", FUN = focal_na_fill_modal)
st_dimensions(obj_focal_filled)
## from to offset delta refsys values x/y
## x 284 306 636540 10 WGS 84 / UTM zone 31N NULL [x]
## y 95 104 5654200 -10 WGS 84 / UTM zone 31N NULL [y]
## t 1 168 NA NA POSIXct 2020-01-06,...,2025-10-10
st_dimensions(obj_classified)
## from to offset delta refsys values x/y
## x 284 306 636540 10 WGS 84 / UTM zone 31N NULL [x]
## y 95 104 5654200 -10 WGS 84 / UTM zone 31N NULL [y]
## t 1 168 NA NA POSIXct 2020-01-06,...,2025-10-10
obj_focal_filled <-
obj_focal_filled |>
st_crop(AOI)
plot(obj_classified)
plot(obj_focal_filled)
# Extract classes into a tibble
obj_focal_filled <- obj_focal_filled - 1 # convert 1-2 to 0-1
classified_df_filled <- as_tibble(obj_focal_filled[AOI])
# Add a month column
classified_df_filled <- classified_df_filled %>%
mutate(month = floor_date(t, "month"))
classified_df_filled_monthly_median <- classified_df_filled %>%
group_by(month, x, y) %>% # group by pixel location + month
summarize(Class = median(X, na.rm = TRUE), .groups = "drop")
# Convert -Inf values to NA
classified_df_filled_monthly_median[classified_df_filled_monthly_median == -Inf] <- NA
# Convert back to stars
classified_filled_monthly_median <- st_as_stars(classified_df_filled_monthly_median, dims = c("x", "y", "month"))
plot(classified_filled_monthly_median)
# Existing data
existing_data <- classified_filled_monthly_median[[1]]
current_months <- st_get_dimension_values(classified_filled_monthly_median, "month")
# Create full monthly sequence
full_months <- seq(from = min(current_months), to = max(current_months), by = "month")
# Prepare new array with NAs for missing months
new_dims <- c(dim(classified_filled_monthly_median)[1:2], length(full_months))
new_array <- array(NA, dim = new_dims)
# Map existing data into correct positions
month_idx <- match(current_months, full_months)
new_array[,,month_idx] <- existing_data
# Start from the original dimensions object
dims <- st_dimensions(classified_filled_monthly_median)
# Update the month dimension's values
dims$month$values <- full_months
# Assign new array to the stars object
classified_filled_monthly_median <- st_as_stars(list(data = new_array))
# Attach updated dimensions
st_dimensions(classified_filled_monthly_median) <- dims
plot(classified_filled_monthly_median)
st_dimensions(classified_filled_monthly_median)
## from to offset delta refsys values x/y
## x 1 23 639370 10 NA NULL [x]
## y 1 10 5653260 -10 NA NULL [y]
## month 1 61 NA NA POSIXct 2020-01-01,...,2025-10-01
# Summarise per month
table_month <-
classified_filled_monthly_median[,,,] |>
as_tibble() |>
group_by(month) |>
summarize(inundation = mean(data, na.rm = TRUE)*100)
table_month_plotting <- table_month[!is.na(table_month$inundation),] # drop NAs for plotting
ggplot(table_month_plotting) + ## check the dates in table to plot a specific single year
aes(x = month, y = inundation) +
ylim(0,100) +
geom_line() + geom_point() +
labs(title = "Inundated area %",
x = "Date",
y = "Inundated area (%)") +
theme_minimal()
# extract time points per month
time_vals1m <- st_get_dimension_values(classified_filled_monthly_median, "month")
time_vals <- 1:dim(classified_filled_monthly_median)[3] # or 1:length(time dimension)
obj_filled_1m <- st_apply(
classified_filled_monthly_median,
MARGIN = c("x", "y"),
FUN = function(ts) {
if (all(is.na(ts))) {
return(ts)
}
approx(
x = seq_along(ts), # numeric indices of time points
y = ts,
xout = seq_along(ts), # interpolate at same time points
method = "linear",
rule = 2
)$y
},
.fname = "month"
)
# fix the broken time dimension in output
obj_filled_1month <- st_set_dimensions(obj_filled_1m, "month", values = time_vals1m)
st_dimensions(classified_filled_monthly_median)
## from to offset delta refsys values x/y
## x 1 23 639370 10 NA NULL [x]
## y 1 10 5653260 -10 NA NULL [y]
## month 1 61 NA NA POSIXct 2020-01-01,...,2025-10-01
st_dimensions(obj_filled_1month)
## from to offset delta refsys values x/y
## month 1 70 NA NA POSIXct 2020-01-01,...,2025-10-01
## x 1 23 639370 10 NA NULL [x]
## y 1 10 5653260 -10 NA NULL [y]
ggplot() +
ggtitle("Per month mean class (Jussila)") +
geom_stars(data = obj_filled_1month["data"]) +
geom_sf(data = AOI, fill = NA, color = "red", linewidth = 1) +
facet_wrap(~month,ncol=12) +
theme_minimal() +
theme(axis.text = element_blank(), # remove axis tick labels
axis.ticks = element_blank(), # remove tick marks
) +
scale_fill_gradientn(
name = "Wetness",
colors = c("tan", "cyan", "blue"),
na.value = "gray",
limits = c(0, 1)
)
# Plot time series per month (gapfilled)
# first summarise for site area:
table_gf1m <- # Gapfilled table per month
obj_filled_1m[,,,] |>
as_tibble() |>
group_by(month) |>
summarize(inundation = mean(data, na.rm = TRUE)*100)
# New column: Flag gapfilled values
table_gf1m$gapfilled <- "monthly median"
table_gf1m$gapfilled[is.na(table_month$inundation)] <- "gapfilled"
table_gf1m$gapfilled <- as.factor(table_gf1m$gapfilled)
# time format to Date for plotting
table_gf1m$date <- as.Date(table_gf1m$month)
ggplot(table_gf1m) + #check the dates in table to plot a single year ([65:74,] -> 2018 in this case)
aes(x = date, y = inundation) + scale_x_date(date_labels = "%b %y") +
geom_line() + geom_point(aes(colour = gapfilled), size=2) +
scale_color_manual(values = c("gapfilled" = "red", "monthly median" = "black"), name=NULL)+
labs(title = "Inundation %",x = "Date",
y = "Inundated area (%)") +
ylim(0,100)+
theme_minimal()
Plotting yearly indicators:
# Using per month gapfilled stars object to calculate statistics
# aggregate stars into yearly layers
## Mean (for each pixel?)
obj_y_mean2 <- aggregate(obj_filled_1month, by = "year",
FUN=mean, na.rm=T)
## Min (for each pixel?)
obj_y_min2 <- aggregate(obj_filled_1month, by = "year",
FUN=min, na.rm=T) # returns Inf for NA areas
obj_y_min2[obj_y_min2 == Inf] <- NA # convert Inf back to NA
## Max (for each pixel?)
obj_y_max2 <- aggregate(obj_filled_1month, by = "year",
FUN=max, na.rm=T) # returns -Inf for NA areas
obj_y_max2[obj_y_max2 == -Inf] <- NA # convert -Inf back to NA
## Plot yearly raster maps
# Inundation (Jussila). Inundated % of time
plot(obj_y_mean2[1], col = viridis(100, option = "D", direction=-1),
breaks = seq(0, 1, length.out = 101), main = "Relative frequency of inundated pixels by Jussila model")
plot(obj_y_mean2[1]*12, col = viridis(100, option = "D", direction=-1),
breaks = seq(0, 12, length.out = 101), main = "Number of inundated months per year by Jussila model")
plot(obj_y_min2[1], col = viridis(100, option = "D", direction=-1),
breaks = seq(0, 1, length.out = 101))
plot(obj_y_max2[1], col = viridis(100, option = "D", direction=-1),
breaks = seq(0, 1, length.out = 101))
# MEAN - summarise values for site area
table_y_mean2 <-
obj_y_mean2[,,,] |>
as_tibble() |>
group_by(time) |>
summarize(inundation = mean(data, na.rm = TRUE)*100)
#summarise min for sites (site mean at minimum wet situation)
table_y_min2 <-
obj_y_min2[,,,] |>
as_tibble() |>
group_by(time) |>
summarize(inundation = mean(data, na.rm = TRUE)*100)
#summarise max for sites (site mean at minimum wet situation)
table_y_max2 <-
obj_y_max2[,,,] |>
as_tibble() |>
group_by(time) |>
summarize(inundation = mean(data, na.rm = TRUE)*100)
# plot percentage of permanently and seasonally wet area.
# Jussila model
ggplot() +
ylim(0,100)+
geom_ribbon(aes(x = table_y_mean2$time, # permanent
ymin = 0,
ymax = table_y_min2$inundation),fill="#1f9e89", alpha= 0.6)+
geom_ribbon(aes(x = table_y_mean2$time, # seasonal
ymin = table_y_min2$inundation,
ymax = table_y_max2$inundation),fill="#6ece58", alpha= 0.3)+ #Temporarily inundated
geom_ribbon(aes(x = table_y_mean2$time, # never
ymin = table_y_max2$inundation,
ymax = 100),fill="#fde725", alpha= 0.2)+ # Never inundated
labs(title = "Permanently and seasonally inundated area (Jussila model)",
x = "Year",
y = "Inundated area %") +
theme_minimal() + theme(legend.position="none") # + scale_x_date(date_breaks = "1 year")
# Extract the layer
x <- obj_y_mean2[1]
# Classify pixel values
vals <- x[[1]]
classified_vals <- ifelse(
is.na(vals), NA,
ifelse(vals == 0, "permanently dry",
ifelse(vals == 1, "permanently inundated",
"seasonally inundated"))
)
# Rebuild a stars object, preserving geometry
classified <- st_as_stars(classified_vals)
st_dimensions(classified) <- st_dimensions(x)
st_crs(classified) <- st_crs(x)
# Convert to factor (only with existing levels)
present_classes <- sort(unique(na.omit(as.vector(classified_vals))))
classified[[1]] <- factor(classified[[1]], levels = present_classes)
# Define full color scheme, but only use what’s needed
full_colors <- c(
"permanently dry" = "tan",
"seasonally inundated" = "lightblue",
"permanently inundated" = "blue4"
)
class_colors <- full_colors[present_classes]
ggplot() +
geom_stars(data = classified, aes(x = x, y = y)) +
scale_fill_manual(values = class_colors, name = "Inundation Class") +
geom_sf(data = AOI, fill = NA, color = "red", size = 0.7) +
coord_sf() +
facet_wrap(~ time, ncol = 3) +
labs(
title = "Inundation Frequency Classification",
x = "Longitude",
y = "Latitude"
) +
theme_minimal() +
theme(
legend.position = "right",
panel.border = element_rect(color = "grey40", fill = NA),
strip.text = element_text(face = "bold")
)